{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "ename": "SyntaxError",
     "evalue": "invalid syntax (<ipython-input-1-e8744d5a7bd4>, line 31)",
     "output_type": "error",
     "traceback": [
      "\u001b[1;36m  File \u001b[1;32m\"<ipython-input-1-e8744d5a7bd4>\"\u001b[1;36m, line \u001b[1;32m31\u001b[0m\n\u001b[1;33m    except Exception,err:\u001b[0m\n\u001b[1;37m                    ^\u001b[0m\n\u001b[1;31mSyntaxError\u001b[0m\u001b[1;31m:\u001b[0m invalid syntax\n"
     ]
    }
   ],
   "source": [
    "# encoding=utf-8\n",
    "import sys\n",
    "import pylab as pl\n",
    "from collections import defaultdict,Counter,OrderedDict\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "# http://www.jb51.net/article/75178.htm\n",
    "from sklearn import datasets\n",
    "\n",
    "# 经纬度计算\n",
    "#计算经纬度坐标距离（米），参考：http://blog.csdn.net/vernice/article/details/46581361\n",
    "from math import radians, cos, sin, asin, sqrt,pi\n",
    "\n",
    "def quickDistance(lon1,lat1,lon2,lat2,max_diff):\n",
    "    # 快速判断两点距离是否超过上限,经纬度差超过max_diff时，直接判定\n",
    "    if abs( lon1 - lon2 )>max_diff or abs( lat1 - lat2 )>max_diff:\n",
    "        return 1\n",
    "    else:\n",
    "        return 0\n",
    "    \n",
    "def calDistance(lon1, lat1, lon2, lat2): \n",
    "    \"\"\" \n",
    "    计算任意两点间的球面距离(m),传参：经度1，纬度1，经度2，纬度2 （十进制度数）  \n",
    "    \"\"\"  \n",
    "    # 转成float类型\n",
    "    try:\n",
    "        lon1 = float(lon1)\n",
    "        lat1 = float(lat1)\n",
    "        lon2 = float(lon2)\n",
    "        lat2 = float(lat2)\n",
    "    except Exception,err:\n",
    "        print >> sys.stderr,'[ERROR] 传参错误！,非数值形式,退出'\n",
    "        return -1\n",
    "    # 空值判断\n",
    "    if lon1 == 0.0 or lat1 == 0.0 or lon2 == 0.0 or lat2 == 0.0:\n",
    "        print >>sys.stderr,'[ERROR] 取值为0，非法，退出'\n",
    "        return -1\n",
    "    # 将十进制度数转化为弧度  \n",
    "    lon1, lat1, lon2, lat2 = map(radians, [lon1, lat1, lon2, lat2])  \n",
    "    # haversine公式  \n",
    "    dlon = lon2 - lon1   \n",
    "    dlat = lat2 - lat1   \n",
    "    a = sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**2  \n",
    "    c = 2 * asin(sqrt(a))   \n",
    "    r = 6371 # 地球平均半径，单位为公里  \n",
    "    return c * r * 1000 \n",
    "\n",
    "# 读取数据 D:\\work\\用户建模画像\\家公司挖掘\\code\\warren.xls\n",
    "print '开始读取数据'\n",
    "# 数据格式：time lon lat adcode geohash\n",
    "#df = pd.read_excel('C:\\Users\\warren\\Desktop\\warren.xlsx')\n",
    "#df = pd.read_excel('C:\\Users\\warren\\Desktop\\warren1.xlsx')\n",
    "df = pd.read_csv('C:\\Users\\warren\\Desktop\\demo.csv')\n",
    "#df = pd.read_csv('C:\\Users\\warren\\Desktop\\\\tmp_location_test1.csv')\n",
    "#df.head()\n",
    "print '提取经纬度信息'\n",
    "# 如果数据量超长，截断\n",
    "max_count = 1000\n",
    "start = 0\n",
    "if len(df)>max_count:\n",
    "    start = max_count\n",
    "#df1= df[df.uid=='BA1C6537-599E-4CBD-8E0C-7E183C4B7450']\n",
    "data = df.iloc[:start,[1,2,3]] # 提取经纬度信息\n",
    "#data = df.iloc[:,[1,2]] # 提取经纬度信息\n",
    "print '数据点信息：编号、lon、lat、geohash\\n',data\n",
    "#exit()\n",
    "points = data.values.tolist()\n",
    "#data = datasets.load_iris()\n",
    "#print data.data\n",
    "#points = [ [i[0],i[1]] for i in data.data[:,:]]\n",
    "#x=raw_input()\n",
    "#points = [[int(eachpoint.split(\"#\")[0]), int(eachpoint.split(\"#\")[1])] for eachpoint in open(\"points\",\"r\")]\n",
    "#==========DBSCAN开始===========\n",
    "# （1）计算每个数据点相邻的数据点\n",
    "#      邻域定义为以该点为中心以边长为2*EPs的网格——右三角矩阵\n",
    "Eps = 200 # 100m的领域\n",
    "#Eps = 0.28\n",
    "#dist_dict = {}# 距离字典,2->3:2.331,key为组合字符串\n",
    "# surroundPoints记录每个点的领域集，对称矩阵\n",
    "surroundPoints = defaultdict(list) # defualt工厂函数，每定义一个key，就会对应生成默认的list类型空值\n",
    "for idx1,point1 in enumerate(points):\n",
    "    for idx2,point2 in enumerate(points):\n",
    "        if (idx1 < idx2):\n",
    "            #tmp_dist = sqrt((point1[0]-point2[0])*(point1[0]-point2[0])+(point1[1]-point2[1])*(point1[1]-point2[1]))\n",
    "            tmp_dist = calDistance(point1[0],point1[1],point2[0],point2[1])\n",
    "            #dist_dict['%s->%s'] = calDistance(point1[0],point1[1],point2[0],point2[1])\n",
    "            #dist_dict['%s->%s'%(idx1,idx2)] = tmp_dist\n",
    "            if(tmp_dist<=Eps):\n",
    "                surroundPoints[idx1].append(idx2)\n",
    "                surroundPoints[idx2].append(idx1)\n",
    "# （2）计算核心点：定义邻域内相邻的数据点的个数大于4的为核心点\n",
    "#MinPts = 5\n",
    "MinPts = 3\n",
    "corePointIdx = [pointIdx for pointIdx,surPointIdxs in surroundPoints.iteritems() if len(surPointIdxs)>=MinPts]\n",
    "# （3）计算边界点：邻域内包含某个核心点的非核心点，定义为边界点\n",
    "borderPointIdx = []\n",
    "for pointIdx,surPointIdxs in surroundPoints.iteritems():\n",
    "    if (pointIdx not in corePointIdx):\n",
    "        for onesurPointIdx in surPointIdxs:\n",
    "            if onesurPointIdx in corePointIdx:\n",
    "                borderPointIdx.append(pointIdx)\n",
    "                break\n",
    "# （4）识别噪音：噪音点既不是边界点也不是核心点\n",
    "noisePointIdx = [pointIdx for pointIdx in range(len(points)) if pointIdx not in corePointIdx and pointIdx not in borderPointIdx]\n",
    "# 核心点、边界点、噪音点\n",
    "#corePoint = [points[pointIdx] for pointIdx in corePointIdx] \n",
    "#borderPoint = [points[pointIdx] for pointIdx in borderPointIdx]\n",
    "#noisePoint = [points[pointIdx] for pointIdx in noisePointIdx]\n",
    "groups = [idx for idx in range(len(points))]\n",
    "# （5）聚类:①先处理核心点②并入边界点\n",
    "# ①各个核心点与其邻域内的所有核心点放在同一个簇中\n",
    "for pointidx,surroundIdxs in surroundPoints.iteritems():\n",
    "    for oneSurroundIdx in surroundIdxs:\n",
    "        if (pointidx in corePointIdx and oneSurroundIdx in corePointIdx and pointidx < oneSurroundIdx):\n",
    "            # 将编号大的节点数值改成同类小编号\n",
    "            for idx in range(len(groups)):\n",
    "                if groups[idx] == groups[oneSurroundIdx]:\n",
    "                    groups[idx] = groups[pointidx]\n",
    "# ②边界点跟其邻域内的某个核心点放在同一个簇中\n",
    "for pointidx,surroundIdxs in surroundPoints.iteritems():\n",
    "    for oneSurroundIdx in surroundIdxs:\n",
    "        if (pointidx in borderPointIdx and oneSurroundIdx in corePointIdx):\n",
    "            # 边界节点的存储的数值改成邻近核心节点编号\n",
    "            groups[pointidx] = groups[oneSurroundIdx]\n",
    "            break\n",
    "# \n",
    "print '邻接矩阵：'\n",
    "print '\\n'.join(['\\t'.join([str(i)]+[str(j) for j in v]) for i,v in surroundPoints.iteritems()])\n",
    "print '核心点集：',corePointIdx\n",
    "print '边界点集：',borderPointIdx\n",
    "print '噪声点集：',noisePointIdx\n",
    "print '聚类结果：',groups\n",
    "#print '\\n'.join(['\\t'.join([[str(i),str(j)] for i,j in enumerate(groups)])])\n",
    "\n",
    "out = [[k,v] for k,v in Counter(groups).iteritems() if v > 1]\n",
    "print out\n",
    "# 计算聚类簇分组\n",
    "result = []\n",
    "for i in range(len(groups)):\n",
    "    if i in noisePointIdx:\n",
    "        # 忽略噪声点\n",
    "        continue\n",
    "    if i == groups[i]:\n",
    "        # 类中第一个节点，格式:[中心点,数目,成员点,密度,平均点,平均距离]\n",
    "        result.append([i,1,[i],-1,[],-1])\n",
    "    elif i > groups[i]:\n",
    "        # 类中后续节点\n",
    "        for j in range(len(result)):\n",
    "            if result[j][0] == groups[i]:\n",
    "                result[j][1] += 1\n",
    "                result[j][2].append(i)\n",
    "                break\n",
    "    else: # i < groups[i]不存在\n",
    "        pass\n",
    "result = filter(lambda x:x[1]>1,result) # 过滤\n",
    "#import copy\n",
    "#result = copy.deepcopy(result_raw)\n",
    "result.sort(key=lambda c: c[1],reverse=True) # 按照类节点多少排序\n",
    "cluster_num = len(result)\n",
    "#print '聚类分组：','\\n'.join(['\\t'.join([repr(j) for j in i]) for i in result])\n",
    "# 计算各类中心点及密度、平均距离\n",
    "for i in range(len(result)):\n",
    "    xy_list = [[j]+points[j] for j in result[i][2]] # [id lon lat]\n",
    "    num = len(xy_list)\n",
    "    x_min = min([j[1] for j in xy_list])\n",
    "    x_max = max([j[1] for j in xy_list])\n",
    "    y_min = min([j[2] for j in xy_list])\n",
    "    y_max = max([j[2] for j in xy_list])\n",
    "    x_avg = sum(([j[1] for j in xy_list]))/num\n",
    "    y_avg = sum(([j[2] for j in xy_list]))/num\n",
    "    # 找代表点（离质心近，且代表大多数点，距离并非首要因素，平均距离最小，消除噪声点的影响）\n",
    "    dist_list = [] # 质心到各点的距离\n",
    "    dist_avg_list = [] # 各点作为中心点时的簇平均距离\n",
    "    dist_min = calDistance(x_min,y_min,x_max,y_max) # 与平均点最小距离初始化为矩形框对角线\n",
    "    dist_i = xy_list[0][0] # 与质心最近的点编号，初始化为第一个点\n",
    "    for j in range(num):\n",
    "        cur_dist = calDistance(x_avg,y_avg,xy_list[j][1],xy_list[j][2])\n",
    "        dist_list.append(cur_dist)\n",
    "        if cur_dist <= dist_min:\n",
    "            dist_min = cur_dist\n",
    "            #dist_i = xy_list[j][0]\n",
    "        # 计算到其它点的平均距离\n",
    "        cur_dist_sum = 0\n",
    "        for k in range(num):\n",
    "            if k == j:\n",
    "                continue\n",
    "            '''\n",
    "            elif k < j:\n",
    "                cur_dist_sum += dist_dict['%s->%s'%(xy_list[k][0],xy_list[j][0])]\n",
    "            else: # k>j\n",
    "                cur_dist_sum += dist_dict['%s->%s'%(xy_list[j][0],xy_list[k][0])]\n",
    "            '''\n",
    "            # 时间换空间,每次都计算距离,不再使用字典\n",
    "            cur_dist_sum += calDistance(xy_list[k][1],xy_list[k][2],xy_list[j][1],xy_list[j][2])\n",
    "        dist_avg_list.append(cur_dist_sum/(num-1))\n",
    "    print >>sys.stdout,'平均距离：%s'%(repr(dist_avg_list))\n",
    "    # 计算簇最小平均距离\n",
    "    dist_i = xy_list[dist_avg_list.index(min(dist_avg_list))][0]\n",
    "    # 计算平均距离\n",
    "    dist_avg = sum(dist_list)/num\n",
    "    # 计算密度(欧氏距离近似)\n",
    "    area = (x_max - x_min)*(y_max - y_min)\n",
    "    if area <= 0:\n",
    "        print >>sys.stderr,'面积计算出错，找不到矩形框,x=[%s,%s],y=[%s,%s]'%(x_min,x_max,y_min,y_max)\n",
    "        phou = -1\n",
    "    else:\n",
    "        phou = result[i][1]/area # 计算簇密度\n",
    "    result[i][0] = dist_i\n",
    "    result[i][3] = phou\n",
    "    result[i][4] = [x_avg,y_avg]\n",
    "    result[i][5] = dist_avg\n",
    "centerPointsIdx = [i[0] for i in result]\n",
    "avgPoints = [i[4] for i in result]\n",
    "print '%s个点包含%s个噪声点,%s个核心点,%s个边界点,最后一共分成%s类'%(len(groups),len(noisePointIdx),len(corePointIdx),len(borderPointIdx),cluster_num)\n",
    "print '聚类分组：','\\n'.join(['\\t'.join([repr(j) for j in i]) for i in result])\n",
    "print '中心点：'\n",
    "print centerPointsIdx\n",
    "print '质心点：'\n",
    "print avgPoints\n",
    "print '#'*100\n",
    "#################组装聚类结果########################\n",
    "# 顶层结构：[uid all cluster_num cluster_info]\n",
    "# cluster_info结构：\\001分隔，每组格式：[id num density center_point centroid_point avg_dist point_list]\n",
    "out_list = ['uid',len(points),cluster_num,'-']\n",
    "cluster_list = []\n",
    "for i in range(cluster_num):\n",
    "    # 每个类格式：[中心点,数目,成员点,密度,平均点,平均距离]\n",
    "    cluster_list.append([str(i+1),str(result[i][1]),'%.3f'%(result[i][3]),'%s,%s,%s'%(points[result[i][0]][0],points[result[i][0]][1],points[result[i][0]][2]),'%s,%s'%(result[i][4][0],result[i][4][1]),'%.3f'%(result[i][5]),'&'.join(['%s,%s,%s'%(points[j][0],points[j][1],points[j][2]) for j in result[i][2]])])\n",
    "out_list[-1] = '\\001'.join(['\\002'.join(i) for i in cluster_list])\n",
    "print out_list\n",
    "#############################可视化部分#####################################\n",
    "# json形式输出，js页面可视化\n",
    "import json\n",
    "tmp_data = data.values.tolist()\n",
    "out = []\n",
    "for i in tmp_data:\n",
    "    out.append(dict(zip(['lng','lat'],i)))\n",
    "out_str = json.dumps(out,indent=4)\n",
    "import re\n",
    "f=open(r\"C:\\Users\\warren\\Desktop\\cluster_update.html\", \"r+\")\n",
    "open('C:\\Users\\warren\\Desktop\\cluster_update1.html', 'w').write(re.sub(r'\\$data', out_str, f.read()))\n",
    "import webbrowser\n",
    "webbrowser.open('C:\\Users\\warren\\Desktop\\cluster_update1.html')\n",
    "# 附：高德地图api开发指南；http://www.cnblogs.com/milkmap/p/3707711.html\n",
    "\n",
    "#\"\"\"\n",
    "# 绘制各个聚类簇,按照形状、颜色依次组合\n",
    "colors = ['r','b','g','w','y','c','m'];shapes = ['o','s','*','d','p','^','v','1','2','3','4']\n",
    "for i in range(cluster_num):\n",
    "    # 中心点编号,数目,成员,密度,平均点\n",
    "    pl.plot([points[j][0] for j in result[i][2]],[points[j][1] for j in result[i][2]],'%s%s'%(shapes[i/len(shapes)],colors[i%len(colors)]))\n",
    "    # [id num density dist]\n",
    "    pl.annotate('%s,%s,%.2f,%.2f'%(i+1,result[i][1],result[i][3],result[i][5]), xy=result[i][4], xytext=result[i][4],arrowprops=dict(facecolor='blue', shrink=0.03))\n",
    "    #pl.annotate('cluster:%s,%s,%s'%(i,result[i][1],result[i][3]), xy=result[i][4], xytext=result[i][4]+0.5,arrowprops=dict(facecolor='blue', shrink=0.03))\n",
    "# 打印中心点,平均点\n",
    "pl.plot([points[i][0] for i in centerPointsIdx], [points[i][1] for i in centerPointsIdx], 'xc',label='0',markersize=20)\n",
    "pl.plot([i[0] for i in avgPoints], [i[1] for i in avgPoints], '+m',label='0',markersize=20)\n",
    "# 打印噪音点，黑色\n",
    "pl.plot([eachpoint[0] for eachpoint in noisePoint], [eachpoint[1] for eachpoint in noisePoint], 'ok')  \n",
    "pl.xlabel('X');pl.ylabel('Y')\n",
    "pl.title('DBSCAN Cluster Result:%s'%(cluster_num))\n",
    "pl.grid()\n",
    "pl.show()\n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "ename": "SyntaxError",
     "evalue": "Missing parentheses in call to 'print'. Did you mean print('开始读取数据')? (<ipython-input-2-4ea83eb852c9>, line 7)",
     "output_type": "error",
     "traceback": [
      "\u001b[1;36m  File \u001b[1;32m\"<ipython-input-2-4ea83eb852c9>\"\u001b[1;36m, line \u001b[1;32m7\u001b[0m\n\u001b[1;33m    print '开始读取数据'\u001b[0m\n\u001b[1;37m                 ^\u001b[0m\n\u001b[1;31mSyntaxError\u001b[0m\u001b[1;31m:\u001b[0m Missing parentheses in call to 'print'. Did you mean print('开始读取数据')?\n"
     ]
    }
   ],
   "source": [
    "import pylab as pl\n",
    "from collections import defaultdict,Counter\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "# 读取数据 D:\\work\\用户建模画像\\家公司挖掘\\code\\warren.xls\n",
    "print '开始读取数据'\n",
    "# 数据格式：time lon lat adcode geohash\n",
    "df = pd.read_excel('C:\\Users\\warren\\Desktop\\warren.xlsx')\n",
    "#df.head()\n",
    "print '提取经纬度信息'\n",
    "data = df.iloc[:,[1,2]] # 提取经纬度信息\n",
    "# json形式输出，js页面可视化\n",
    "import json\n",
    "tmp_data = data.values.tolist()\n",
    "out = []\n",
    "for i in tmp_data:\n",
    "    out.append(dict(zip(['lng','lat'],i)))\n",
    "out_str = json.dumps(out,indent=4)\n",
    "import re\n",
    "f=open(r\"C:\\Users\\warren\\Desktop\\cluster_update.html\", \"r+\")\n",
    "open('C:\\Users\\warren\\Desktop\\cluster_update1.html', 'w').write(re.sub(r'\\$data', out_str, f.read()))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "ename": "SyntaxError",
     "evalue": "invalid syntax (<ipython-input-3-6ab869a8b2eb>, line 7)",
     "output_type": "error",
     "traceback": [
      "\u001b[1;36m  File \u001b[1;32m\"<ipython-input-3-6ab869a8b2eb>\"\u001b[1;36m, line \u001b[1;32m7\u001b[0m\n\u001b[1;33m    print locals()\u001b[0m\n\u001b[1;37m               ^\u001b[0m\n\u001b[1;31mSyntaxError\u001b[0m\u001b[1;31m:\u001b[0m invalid syntax\n"
     ]
    }
   ],
   "source": [
    "import sys\n",
    "import math\n",
    "\n",
    "# 经纬度计算\n",
    "def calculateDistance(lon1,lat1,lon2,lat2):\n",
    "    # 计算两个坐标点之间的距离 (lon1,lat1),(lon2,lat2)\n",
    "    print locals()\n",
    "    _earth_r = 6371.393 # 地球半径\n",
    "    _pai = 3.1415926 # π值\n",
    "    # 类型转换\n",
    "    for var in ('lon1','lat1','lon2','lat2'):\n",
    "        # 类型转换\n",
    "        if type(eval(var)) != type(0.0):\n",
    "            try:\n",
    "                tmp = eval('float(var)')\n",
    "                eval(var+'='+tmp)\n",
    "                #eval(var + '='+ 'float(' + var + ')')\n",
    "                #eval('%s=float(%s)'%(var))\n",
    "                #eval(var) = float(eval(var))\n",
    "                #eval('var = double(var)',locals())\n",
    "            except Exception,err:\n",
    "                print >> sys.stderr,'[ERROR] 传参错误！(%s),%s非数值形式,退出'%(err,var)\n",
    "                return 0\n",
    "        # 空值判断\n",
    "        if eval(var) == 0.0:\n",
    "            print >>sys.stderr,'[ERROR] %s取值为0，非法，退出'%(var)\n",
    "            return 0\n",
    "    # 计算球面距离\n",
    "    #dist = _earth_r*Math.acos(Math.sin(x1*pai/180)*Math.sin(x2*pai/180)+Math.cos(x1*pai/180)*Math.cos(x2*pai/180)*Math.cos(y1*pai/180-y2*pai/180));\n",
    "    dist = _earth_r*math.acos(math.sin(lat1*_pai/180)*math.sin(lat2*_pai/180)+math.cos(lat1*_pai/180)*math.cos(lat2*_pai/180)*math.cos(lon1*_pai/180-lon2*_pai/180))\n",
    "    return dist\n",
    "\n",
    "if __name__ == '__main__':\n",
    "    calculateDistance(\"30.579249\", \"103.959226\", \"tes0.578909\", \"103.958961\") # 非法取值\n",
    "    calculateDistance(\"30.579249\", \"103.959226\", \"0.0\", \"103.958961\") # 一个空值\n",
    "    calculateDistance(\"30.579249\", \"103.959226\", \"30.579249\", \"103.959226\") # 两个坐标相同\n",
    "    calculateDistance(\"30.579249\", \"103.959226\", \"30.578909\", \"103.958961\")\n",
    "    calculateDistance(\"30.579249\", \"103.959226\", \"30.636895\", \"104.029880\")\n",
    "    calculateDistance(\"40.0\", \"116.5\", \"31.0\", \"121.5\")\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 197,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-1\n",
      "-1\n",
      "9.50416922569e-05\n",
      "0.0309671075332\n",
      "8.03913647571\n",
      "738.614952468\n",
      "896.532784703\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[ERROR] 传参错误！,非数值形式,退出\n",
      "[ERROR] 取值为0，非法，退出\n"
     ]
    }
   ],
   "source": [
    "import sys\n",
    "from math import *\n",
    "\n",
    "# 经纬度计算——注：有问题，个别ｃａｓｅ分母为０，无法运行\n",
    "def calDistance(lon1,lat1,lon2,lat2):\n",
    "    # 计算两个坐标点之间的距离 (lon1,lat1),(lon2,lat2)\n",
    "    # 代码参考：http://www.cnblogs.com/yejingcn/p/4537863.html?utm_source=tuicool&utm_medium=referral\n",
    "    #print locals()\n",
    "    # 类型转换\n",
    "    try:\n",
    "        lon1 = float(lon1)\n",
    "        lat1 = float(lat1)\n",
    "        lon2 = float(lon2)\n",
    "        lat2 = float(lat2)\n",
    "    except Exception,err:\n",
    "        print >> sys.stderr,'[ERROR] 传参错误！,非数值形式,退出'\n",
    "        return -1\n",
    "    # 空值判断\n",
    "    if lon1 == 0.0 or lat1 == 0.0 or lon2 == 0.0 or lat2 == 0.0:\n",
    "        print >>sys.stderr,'[ERROR] 取值为0，非法，退出'\n",
    "        return -1\n",
    "    # 计算球面距离\n",
    "    #dist = _earth_r*Math.acos(Math.sin(x1*pai/180)*Math.sin(x2*pai/180)+Math.cos(x1*pai/180)*Math.cos(x2*pai/180)*Math.cos(y1*pai/180-y2*pai/180));\n",
    "    #dist = _earth_r*math.acos(math.sin(lat1*_pai/180)*math.sin(lat2*_pai/180)+math.cos(lat1*_pai/180)*math.cos(lat2*_pai/180)*math.cos(lon1*_pai/180-lon2*_pai/180))\n",
    "    ra = 6378.140  # 赤道半径 (km)\n",
    "    rb = 6356.755  # 极半径 (km)\n",
    "    flatten = (ra - rb) / ra  # 地球扁率\n",
    "    rad_lat_A = radians(lat1)\n",
    "    rad_lng_A = radians(lon1)\n",
    "    rad_lat_B = radians(lat2)\n",
    "    rad_lng_B = radians(lon2)\n",
    "    pA = atan(rb / ra * tan(rad_lat_A))\n",
    "    pB = atan(rb / ra * tan(rad_lat_B))\n",
    "    xx = acos(sin(pA) * sin(pB) + cos(pA) * cos(pB) * cos(rad_lng_A - rad_lng_B))\n",
    "    c1 = (sin(xx) - xx) * (sin(pA) + sin(pB)) ** 2 / cos(xx / 2) ** 2\n",
    "    c2 = (sin(xx) + xx) * (sin(pA) - sin(pB)) ** 2 / sin(xx / 2) ** 2\n",
    "    dr = flatten / 8 * (c1 - c2)\n",
    "    distance = ra * (xx + dr)\n",
    "    return distance\n",
    "\n",
    "if __name__ == '__main__':\n",
    "    print calDistance(\"30.579249\", \"103.959226\", \"tes0.578909\", \"103.958961\") # 非法取值\n",
    "    print calDistance(\"30.579249\", \"103.959226\", \"0.0\", \"103.958961\") # 一个空值\n",
    "    print calDistance(\"30.579249\", \"103.959226\", \"30.579249\", \"103.959226\") # 两个坐标相同\n",
    "    print calDistance(\"30.579249\", \"103.959226\", \"30.578909\", \"103.958961\")\n",
    "    print calDistance(\"30.579249\", \"103.959226\", \"30.636895\", \"104.029880\")\n",
    "    print calDistance(\"40.0\", \"116.5\", \"31.0\", \"121.5\")\n",
    "    print calDistance('118.796877','32.060255','116.407395','39.904211') # 南京到北京的距离896.533 km\n",
    "    \n",
    "\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-1\n",
      "-1\n",
      "0.0\n",
      "30.8457080039\n",
      "8007.82935937\n",
      "736563.355211\n",
      "898211.728891\n",
      "200.15086796\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[ERROR] 传参错误！,非数值形式,退出\n",
      "[ERROR] 取值为0，非法，退出\n"
     ]
    }
   ],
   "source": [
    "#计算经纬度坐标距离，参考：http://blog.csdn.net/vernice/article/details/46581361\n",
    "import sys\n",
    "from math import radians, cos, sin, asin, sqrt ,pi \n",
    "#import math\n",
    "  \n",
    "def calDistance(lon1, lat1, lon2, lat2): # 经度1，纬度1，经度2，纬度2 （十进制度数）  \n",
    "    \"\"\" \n",
    "    Calculate the great circle distance between two points  \n",
    "    on the earth (specified in decimal degrees) \n",
    "    \"\"\"  \n",
    "    # 转成float类型\n",
    "    try:\n",
    "        lon1 = float(lon1)\n",
    "        lat1 = float(lat1)\n",
    "        lon2 = float(lon2)\n",
    "        lat2 = float(lat2)\n",
    "    except Exception,err:\n",
    "        print >> sys.stderr,'[ERROR] 传参错误！,非数值形式,退出'\n",
    "        return -1\n",
    "    # 空值判断\n",
    "    if lon1 == 0.0 or lat1 == 0.0 or lon2 == 0.0 or lat2 == 0.0:\n",
    "        print >>sys.stderr,'[ERROR] 取值为0，非法，退出'\n",
    "        return -1\n",
    "    # 将十进制度数转化为弧度  \n",
    "    lon1, lat1, lon2, lat2 = map(radians, [lon1, lat1, lon2, lat2])  \n",
    "    # haversine公式  \n",
    "    dlon = lon2 - lon1   \n",
    "    dlat = lat2 - lat1   \n",
    "    a = sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**2  \n",
    "    c = 2 * asin(sqrt(a))   \n",
    "    r = 6371 # 地球平均半径，单位为公里  \n",
    "    return c * r * 1000 \n",
    "\n",
    "if __name__ == '__main__':\n",
    "    print calDistance(\"30.579249\", \"103.959226\", \"tes0.578909\", \"103.958961\") # 非法取值\n",
    "    print calDistance(\"30.579249\", \"103.959226\", \"0.0\", \"103.958961\") # 一个空值\n",
    "    print calDistance(\"30.579249\", \"103.959226\", \"30.579249\", \"103.959226\") # 两个坐标相同\n",
    "    print calDistance(\"30.579249\", \"103.959226\", \"30.578909\", \"103.958961\")\n",
    "    print calDistance(\"30.579249\", \"103.959226\", \"30.636895\", \"104.029880\")\n",
    "    print calDistance(\"40.0\", \"116.5\", \"31.0\", \"121.5\")\n",
    "    print calDistance('118.796877','32.060255','116.407395','39.904211') # 南京到北京的距离896.533 km\n",
    "    R = 6371000\n",
    "    print 0.0018*pi*R/180"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 244,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0, 1, [0]], [1, 8, [1, 2, 3, 6, 10, 12, 14, 19]], [4, 1, [4]], [5, 1, [5]], [7, 3, [7, 15, 17]], [8, 1, [8]], [9, 1, [9]], [11, 1, [11]], [13, 1, [13]], [16, 1, [16]], [18, 1, [18]], [20, 1, [20]]]\n",
      "[[1, 8, [1, 2, 3, 6, 10, 12, 14, 19]], [7, 3, [7, 15, 17]]]\n"
     ]
    }
   ],
   "source": [
    "groups = [0, 1, 1, 1, 4, 5, 1, 7, 8, 9, 1, 11, 1, 13, 1, 7, 16, 7, 18, 1, 20]\n",
    "#print noisePointIdx\n",
    "# 计算聚类簇分组\n",
    "result = []\n",
    "for i in range(len(groups)):\n",
    "    #if i in noisePointIdx:\n",
    "    #    # 忽略噪声点\n",
    "    #    continue\n",
    "    if i == groups[i]:\n",
    "        # 类中第一个节点，格式:[中心点,数目,成员点,密度,平均点,平均距离]\n",
    "        result.append([i,1,[i]])\n",
    "    elif i > groups[i]:\n",
    "        # 类中后续节点\n",
    "        for j in range(len(result)):\n",
    "            if result[j][0] == groups[i]:\n",
    "                result[j][1] += 1\n",
    "                result[j][2].append(i)\n",
    "                break\n",
    "    else: # i < groups[i]不存在\n",
    "        pass\n",
    "print result\n",
    "result= filter(lambda x:x[1]>1,result)\n",
    "print result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
